###Democratization boost or bust? Electoral turnout after democratic transitions###
###Created by Roman-Gabriel Olar on 01/05/2023###
###Last modified on 03/07/2023###


all_pack<-c('skimr', 'fabricatr', 'tidyverse', 'haven',  'xlsx',  'pastecs', 'mgcv', 
            'itsadug', 'lme4', 'broom', 'sjPlot', 'gridExtra', 
            'grid', 'truncnorm', 'lattice', 'merTools', 'rlist', 'texreg', 'ggrepel')
lapply(all_pack, library, character.only=TRUE)
library(texreg)
library(haven)          
library(BaM)
library(pacman)
library(mgcv)
library(ggplot2)
library(dplyr)
library(mgcv)
library(tidymv)
library(foreign)
library(ggthemes)
library(officer)
library(broom)
library(stargazer)
setwd("Put your own folder directory")

df<-read_dta("gam_micro.dta")
df$wave=as.factor(df$wave)
df$ccode=as.factor(df$ccode)
df$educ_combi=as.factor(df$educ_combi)
df$year2=as.factor(df$year)
df$data2=as.factor(df$data)
df$ccode_year=as.factor(df$ccode_year)
df$coh_tab=as.factor(df$coh_tab)
df$elec_type=as.factor(df$elec_type)
df$idea=as.factor(df$idea)
df$conc <- ifelse(is.na(df$conc), 0,df$conc)
df$conc=as.factor(df$conc)
df$e_regionpol=as.factor(df$e_regionpol)


###Model A3.5 from Table A3.3###
gam_crt <- bam(turnout ~ s(age_yr_crt) + age + age2 + 
             educ_combi + working + female + l_wbpopln + l_ln_gdppc + proportional + comp_vote +
             majoritarian + conc +  data2 + year2 + e_regionpol + s(ccode, bs = 're'), family = binomial, data=df)
summary(gam_crt)
set_theme(base = theme_light())

###Figure 2 from main text###
plot(gam_crt, select = 1, xlab=" Age at civil resistance transition", ylab="Probability to turnout to vote") 
abline(h=0, col="red", lty = 2) 
axis(1, at = seq(0, 80, by = 10), labels = c("0", "10", "20", "30", "40", "50", "60", "70","80" ))
axis(2, at = seq(-.9, .4, by = .1), labels = c("-0.9","-0.8", "-0.7", "-0.6", "-0.5", "-0.4", "-0.3", "-0.2", "-0.1", "", "0.1", "0.2", "0.3","0.4"))

###Model A3.6 from Table A3.3###
gam_gwf <- bam(turnout ~ s(age_yr_gwf) + age + age2 + 
              educ_combi + working + female + l_wbpopln + l_ln_gdppc + proportional + comp_vote +
              majoritarian + conc +  data2 + year2 + e_regionpol + s(ccode, bs = 're'), family = binomial, data=df)
summary(gam_gwf)

###Figure 3 from main text###
plot(gam_gwf, select = 1, xlab=" Age at other transition", ylab="Probability to turnout to vote") 
abline(h=0, col="red", lty = 2) 
axis(1, at = seq(0, 80, by = 10), labels = c("0", "10", "20", "30", "40", "50", "60", "70","80" ))
axis(2, at = seq(-.4, .5, by = .1), labels = c("-0.4", "-0.3", "-0.2", "-0.1", "", "0.1", "0.2","0.3","0.4","0.5"))


###Model A3.7 from Table A3.3###

###Recode as 0 those from opposing transition type country
df$age_yr_crt <- ifelse(is.na(df$age_yr_crt), 0,df$age_yr_crt)
df$age_yr_gwf <- ifelse(is.na(df$age_yr_gwf), 0,df$age_yr_gwf)

gam_newd <- bam(turnout ~ s(age_yr_crt) + s(age_yr_gwf) + age + age2 + 
             educ_combi + working + female + l_wbpopln + l_ln_gdppc + proportional + comp_vote +
             majoritarian + conc +  data2 + year2 + e_regionpol + s(ccode, bs = 're'), family = binomial, data=df)
summary(gam_newd)

###Figure A3.1###
plot(gam_newd, select = 1, xlab=" Age at civil resistance transition", ylab="Probability to turnout to vote") 
abline(h=0, col="red", lty = 2) 
axis(1, at = seq(0, 80, by = 10), labels = c("0", "10", "20", "30", "40", "50", "60", "70","80" ))
axis(2, at = seq(-.7, .3, by = .1), labels = c("-0.7","-0.6", "-0.5","-0.4", "-0.3", "-0.2", "-0.1", "", "0.1", "0.2","0.3"))

###Figure A3.2###
plot(gam_newd, select = 2, xlab=" Age at other transition", ylab="Probability to turnout to vote") 
abline(h=0, col="red", lty = 2) 
axis(1, at = seq(0, 80, by = 10), labels = c("0", "10", "20", "30", "40", "50", "60", "70","80" ))
axis(2, at = seq(-.8, .3, by = .1), labels = c("-0.8", "-0.7","-0.6", "-0.5","-0.4", "-0.3", "-0.2", "-0.1", "", "0.1", "0.2","0.3"))

###Model citizens from established democracies

df2<-read_dta("Micro_Analysis.dta")
df2$wave=as.factor(df2$wave)
df2$ccode=as.factor(df2$ccode)
df2$educ_combi=as.factor(df2$educ_combi)
df2$year2=as.factor(df2$year)
df2$data2=as.factor(df2$data)
df2$ccode_year=as.factor(df2$ccode_year)
df2$coh_tab=as.factor(df2$coh_tab)
df2$elec_type=as.factor(df2$elec_type)
df2$idea=as.factor(df2$idea)
df2$conc <- ifelse(is.na(df2$conc), 0,df2$conc)
df2$conc=as.factor(df2$conc)
df2$e_regionpol=as.factor(df2$e_regionpol)

###Recode data (new=0 if age<0 and established democracies=0)
df2$age_yr_crt <- ifelse(is.na(df2$age_yr_crt), 0,df2$age_yr_crt)
df2$age_yr_gwf <- ifelse(is.na(df2$age_yr_gwf), 0,df2$age_yr_gwf)
df2$age_yr_crt <- ifelse(df2$age_yr_crt < 0, 0, df2$age_yr_crt)
df2$age_yr_gwf <- ifelse(df2$age_yr_gwf < 0, 0, df2$age_yr_gwf)
df2$age_yr_crt <- ifelse(df2$age_yr_crt > 80, 80, df2$age_yr_crt)
df2$age_yr_gwf <- ifelse(df2$age_yr_gwf > 80, 80, df2$age_yr_gwf)

###Model A3.8###

gam_all <-bam(turnout ~ s(age_yr_crt) + s(age_yr_gwf) + age + age2 + 
           educ_combi + working + female + l_wbpopln + l_ln_gdppc + proportional + comp_vote +
           majoritarian + conc +  data2 + year2 + e_regionpol + s(ccode, bs = 're'), family = binomial, data=df2)

###Figure A3.3###
plot(gam_all, select = 1, xlab=" Age at civil resistance transition", ylab="Probability to turnout to vote") 
abline(h=0, col="red", lty = 2) 
axis(1, at = seq(0, 80, by = 10), labels = c("0", "10", "20", "30", "40", "50", "60", "70","80" ))
axis(2, at = seq(-.8, .3, by = .1), labels = c("-0.8", "-0.7","-0.6", "-0.5","-0.4", "-0.3", "-0.2", "-0.1", "", "0.1", "0.2", "0.3"))

###Figure A3.4###
plot(gam_all, select = 2, xlab=" Age at other transition", ylab="Probability to turnout to vote") 
abline(h=0, col="red", lty = 2) 
axis(1, at = seq(0, 80, by = 10), labels = c("0", "10", "20", "30", "40", "50", "60", "70","80" ))
axis(2, at = seq(-.8, .2, by = .1), labels = c("-0.8", "-0.7","-0.6", "-0.5","-0.4", "-0.3", "-0.2", "-0.1", "", "0.1", "0.2"))

###EXPORT ALL TABLES####

wordreg(list(gam_crt, gam_gwf, gam_newd, gam_all), file = "gam_models.doc")